/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAFNN
  \brief This class represents nucleon-nucleon scattering amplitude, extracted
  from fitting of pp, pn and nn scattering data over a wide range of energies.
  fNN(q)=kNN/4pi*sigmaNN(i+alphaNN)*exp(-betaNN*q^2)
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/08
  \date Last modified: 2020/09/09 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include "TAFNN.h"
#include "TAMath.h"
#include "TAException.h"
#include "TAInterpolate.h"

#define polint TAInterpolateD::PolyInter
#define sqr TAMath::sqr

static const double u = TAMath::uMeV(); // atomic mass unit in MeV
// Nof points: polynomial interpolation order = NP - 1
static const int NP = 6; // NOTE that MOMDIS (fortran) used 6 points uniformly
static const int NPHF = NP/2;

TAFNN::ParOption TAFNN::kParOpt = TAFNN::kHoriuchi;

TAFNN::TAFNN(){}
TAFNN::~TAFNN(){}

void TAFNN::SetParOption(int opt){
  switch(opt){
    case 0: kParOpt = kHoriuchi; break;
    case 1: kParOpt = kLenziRay; break;
    default:
      TAException::Error("TAFNN",
        "SetParOption: illegal input (0: Horiuchi; 1: LenziRay).");
      break;
  } // end switch
} // end member function SetParOption

// P: projectile, T: target
double TAFNN::IsospinAverage(double vpp, double vpn, int zP, int nP, int zT, int nT){
  const double aP = zP + nP, aT = zT + nT;
  if(zP < 0 || zT < 0 || nP < 0 || nT < 0 || aP <= 0 || aT <= 0)
    TAException::Error("TAFNN", "IsospinAverage: One of the following occurs:\
 nzP < 0 || zT < 0 || nP < 0 || nT < 0 || aP <= 0 || aT <= 0");

  return (vpp*(nP*nT+zP*zT) + vpn*(nP*zT+zP*nT)) / (aP*aT);
} // end member function IsospinAverage

/// isospin-averaged nucleon-nucleon scattering cross section
/// for lab energy pEk.
/// \param isPauli 0: free nn scattering; 1: Pauli-blocking corrected
/// \retval in fm^2
double TAFNN::GetSigmaNN(double pEk, int zP, int nP, int zT, int nT,
    bool isPauli){
  double sigmaPP = 0., sigmaPN = 0.; // pp and NP scattering c.s.
  const double b = TAMath::EkPerUToBeta(pEk); // beta=v/c
  const double b2 = b*b, b4 = b2*b2;

  if(!isPauli){
    if(pEk > 10.){
      sigmaPN = -70.67 - 18.18/b + 25.26/b2 + 113.85*b;
      sigmaPP =  13.73 - 15.04/b +  8.76/b2 +  68.67*b4;
    } // end if(pEk > 10.)
    else{
      sigmaPN = (2.73/(0.35*pEk+sqr(1-0.0553*pEk)) +
        17.63/(6.8*pEk+sqr(1+0.344*pEk)))*1.e3;
      sigmaPP =  13.73 - 15.04/b +  8.76/b2 +  68.67*b4;
    } // end else
  } // end if(isPauli)
  else{ // take in Pauli-blocking
    static const double re = 0.17; // effective rho for Pauli-blocking in fm
    sigmaPP =  13.73 - 15.04/b +  8.76/b2 +  68.67*b4;
    sigmaPP *= (1.+7.772*pow(pEk, 0.06)*pow(re, 1.48)) /
      (1.+18.01*pow(re, 1.46));
    sigmaPN = -70.67 - 18.18/b + 25.26/b2 + 113.85*b;
    sigmaPN *= (1.+20.88*pow(pEk, 0.04)*pow(re, 2.02)) /
      (1.+35.86*pow(re, 1.9));
  } /// end the final else
  return IsospinAverage(sigmaPP, sigmaPN, zP, nP, zT, nT) * 0.1; // 0.1: mb->fm^2
} // end member function GetSigmaNN

// various parameterizations exist for alphaNN and betaNN
/// fNN(q)=kNN/4pi*sigmaNN(i+alphaNN)*exp(-betaNN*q^2)
/// alphaNN at arbitrary Ek is interpolated from array alphaNN
double TAFNN::GetAlphaNN(double pEk, int zP, int nP, int zT, int nT){
  if(kHoriuchi == kParOpt) return AlphaBetaHoriuchi(pEk, true);
  return AlphaBetaLenziRay(pEk, zP, nP, zT, nT, true);
} // end member function GetAlphaNN
double TAFNN::GetBetaNN(double pEk, int zP, int nP, int zT, int nT){
  if(kHoriuchi == kParOpt) return AlphaBetaHoriuchi(pEk, false);
  return AlphaBetaLenziRay(pEk, zP, nP, zT, nT, false);
} // end member function GetBetaNN

/// Horiuchi, et al., PRC.75.044607. E: [30, 1000], N=16
double TAFNN::AlphaBetaHoriuchi(double pEk, bool isAlpha){
  if(pEk < 30. || pEk > 1000.) TAException::Error("TAFNN",
    "GetAlphaNN: Ek'd be within [30, 1000] MeV/nucleon for Horiuchi Parameters.");

  static const int N = 16;
  static const double e[N] = { // Elab, in MeV/nucleon
    30., 38., 40., 49., 85., 94., 100., 120., 150., 200.,
    325., 425., 550., 650., 800., 1000.
  };
  static const double a[N] = { // alphaNN, already isospin averaged
    0.87, 0.89, 0.9, 0.94, 1.37, 1.409, 1.435, 1.359, 1.245, 0.953,
    0.305, 0.36, 0.04, -0.095, -0.07, -0.275
  };
  static const double b[N] = { // betaNN, already isospin averaged
    0.685, 0.521, 0.486, 0.39, 0.349, 0.327, 0.322, 0.255, 0.195, 0.131,
    0.075, 0.078, 0.125, 0.16, 0.21, 0.21
  };

  // select a segment to wrap pEk in
  if(isAlpha) return polint(e, a, N, NP, pEk);
  return polint(e, b, N, NP, pEk) / 2.;
} // end member function AlphaBetaHoriuchi

/// Lenzi, et al., PRC.40.2114. E: [10, 94], N=7
/// Ray, et al., PRC.20.1857. E: [100, 2200], N=10
double TAFNN::AlphaBetaLenziRay(double pEk, int zP, int nP, int zT, int nT,
    bool isAlpha){
  if(pEk < 10. || pEk > 2200.) TAException::Error("TAFNN",
    "GetAlphaNN: Ek'd be within [10, 2200] MeV/nucleon for Lenzi & Ray Parameters.");

  static const int N = 17, OFST = 7; // 7 + 10
  static const double e[N] = {
    10., 30., 38., 40., 49., 85., 94.,
    100., 150., 200., 325., 425., 550., 650., 800., 1000., 2200.
  };
  static const double app[N] = {
    0.8, 0.87, 0.89, 0.9, 0.94, 1.0, 1.07,
    1.87, 1.53, 1.15, 0.45, 0.47, 0.32, 0.16, 0.06, -0.09, -0.17
  };
  static const double apn[N] = {
    0.8, 0.87, 0.89, 0.9, 0.94, 1.0, 1.07,
    1.00, 0.96, 0.71, 0.16, 0.25, -0.24, -0.35, -0.20, -0.46, -0.50
  };
  static const double bpp[N] = {
    0., 0., 0., 0., 0., 0., 0.,
    0.66, 0.57, 0.56, 0.26, 0.21, 0.04, 0.07, 0.09, 0.09, 0.12
    // 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25
  };
  static const double bpn[N] = {
    0., 0., 0., 0., 0., 0., 0.,
    0.36, 0.58, 0.68, 0.36, 0.27, 0.085, 0.09, 0.12, 0.12, 0.14
    // 0.25, 0.25, 0.25, 0.25, 0.25, 0.25,  0.25, 0.25, 0.25, 0.25
  };

  if(pEk <= e[6]){ // Lenzi
    if(isAlpha) return IsospinAverage(polint(e, app, OFST, NP, pEk),
      polint(e, apn, OFST, NP, pEk), zP, nP, zT, nT);
    return IsospinAverage(polint(e, bpp, OFST, NP, pEk),
      polint(e, bpn, OFST, NP, pEk), zP, nP, zT, nT);
  } // end if
  else{ // Ray
    // cout << "zP: " << zP << " nP: " << nP << " zT: " << zT << " nT: " << nT << endl;
    // cout << "app: " << polint(e+OFST, app+OFST, N-OFST, NP, pEk) << endl;
    // cout << "apn: " << polint(e+OFST, apn+OFST, N-OFST, NP, pEk) << endl;
    // cout << "bpp: " << polint(e+OFST, bpp+OFST, N-OFST, NP, pEk) << endl;
    // cout << "bpn: " << polint(e+OFST, bpn+OFST, N-OFST, NP, pEk) << endl;
    // getchar();
    if(isAlpha) return IsospinAverage(polint(e+OFST, app+OFST, N-OFST, NP, pEk),
      polint(e+OFST, apn+OFST, N-OFST, NP, pEk), zP, nP, zT, nT);
    return IsospinAverage(polint(e+OFST, bpp+OFST, N-OFST, NP, pEk),
      polint(e+OFST, bpn+OFST, N-OFST, NP, pEk), zP, nP, zT, nT);
  } // end else
} // end member function AlphaBetaLenziRay

#ifdef CATCH_CNOK
TEST_CASE("template class TAInterpolate", "[polint]"){
  double x[4] = {1., 2., 3., 4.}, y[4] = {1., 4., 9., 16.}, dy = 0.;
  // double yy = polint(x, y, 4, -2.5, &dy);
  double yy = TAInterpolateD::PolyInterLag(x, y, 4, -2.5, &dy);
  // TAException::Info("TAInterpolate", "PolyInter: dy is %f", dy);
  CHECK(yy == 6.25);
  CHECK(dy == 0.);
} // end TEST_CASE
#endif
